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Abstract 

A free-energy functional for a crystal that contains both the symmetry conserved and symmetry 
broken parts of the direct pair correlation function is developed. The free-energy functional is 
used to investigate the crystallization of fluids interacting via the inverse power potential ; u{r) = 
e{a/ry^. In agreement with simulation results we find that for n = 12 the freezing is into close 
packed f.c.c structure while for soft repulsions (n < 6) b.c.c phase is more stable. 



PACS numbers: 64.70.dg, 64.70.dm, 05.70.Fh 



When a fluid freezes into a crystalline solid its continuous symmetry of translation and 
rotation is broken into one of the symmetry groups of the Bravais lattices. A crystalline 
solid has a discrete set of vectors Ri such that any function of position, such as one particle 
density satisfies p(r) = p(r+Ri) for all Ri [1]. This set of vectors necessarily forms a Bravais 
lattice. While many metallic systems freeze into b.c.c structure, simple fluids like Ar freeze 
into f.c.c structure [2]. Model fluids interacting via inverse power potentials u{r) = e{a/r)"' 
where e, a and n are potential parameters and r molecular separation, show interesting 
freezing transitions [3-5]. The more repulsive (n > 7) systems freeze into f.c.c structure 
while the soft repulsions n <6 freeze into b.c.c phase. A correct description of the stability 
of the two cubic structures is a stringent test of any statistical mechanical theory based on 
first principle, since atomic arrangements in the two are very different. 

Since 1979 [6] the density functional theory (DFT) has been apphed to the problem of 
crystaUization of a wide variety of fluids [7-8] . Despite its many successes notably with hard 
sphere system, a full theory applicable to all intermolecular potentials has remained elusive 
[7-11]. The direct pair correlation function (DPCF) that appears in the expression of the 
free-energy functional corresponds to the ordered phase and is functional of p(r). When this 
DPCF is replaced by that of the co-existing isotropic liquid [6] or by that of an "effective 
fluid" [12], the free-energy functional becomes approximate and fails to provide a correct 
description of the freezing transition. Attempts to include a term involving three-body 
direct correlation function of the isotropic phase in the free-energy functional have failed to 
improve the situation [9,10]. 

Due to breaking of symmetry at the fluid-solid transition a qualitatively new contribution 
[13,14] to the pair correlation function of a crystal arises . The pair correlation function 
of a crystal has therefore two different contributions; one that maintains the continuous 
symmetry of the Hamiltonian and passes smoothly through the transition and the other 
that breaks it and vanishes at the melting point. In this Letter we develop a method 
to estimate the DPCF of a crystal and use it to construct a free energy functional by 
performing functional integrations in the appropriate domains of density space. We use 
this free-energy functional to investigate the crystallization of fluids interacting via the 
inverse power potential. Potentials of this class have a simple scaling property according 
to which the reduced excess thermodynamic properties depend on a single variable which 
is defined as 7 = (p(T^)(e/A;bT)(^/") = where kj, is the Boltzmann constant and 
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T, temperature. The equation of state and the melting curve of these potentials have been 
extensively investigated by Monte Carlo (MC) simulations for several values of parameter n 
so that "exact" results are available [3-5] for comparison. 

The reduced free-energy A[p\ of an inhomogeneous system is functional of p(r) and is 
written as A[p\ — Aid[p\ + Aex[p\ where Aid — I drp{r)[\n{p{r)A) — 1] is the ideal gas part. 
Here A is the cube of thermal wavelength associated with a particle. The excess part Aex[p\ 
is related with the DPCF of the system as [14] 

where c^°^ is symmetry conserving and c^''^ symmetry broken parts of the DPCF. While 
c^°^ depends on the magnitude of the interparticle separation and is function of the average 
number density po, c^^^ is invariant only under discrete set of translations and is functional 
of p(r). If one chooses a center of mass variable Vf. = (ri + r2)/2 and difference variable 
r = ri — r2, then c^''^ can be written as [15]; 



c 



ri,r2;[p]) = Eexpi(G-re)c(«)(r;[p]) (2) 

G 

where G are reciprocal lattice vectors (R. L. V.). Since the DPCF is real and symmetric 
with respect to interchange of ri and r2, c^^\r) = c^~^\r) and c*^'^^(r) = c^^\—r) 

The Ornstein-Zernike equation and the Roger- Young closure relation [16] have been used 
to calculate c*^°^ (r) and '^'^^q^^^ as a function of intermolecular separation for densities ranging 
from zero to a density above the melting point at small intervals for n=12,6 and 4. This 
method gives values of pair correlation functions which are thermodynamically self consis- 
tent. For c(^)(ri,r2,[p]) one can either generalize the integral equation method [14] or use a 
functional Taylor expansion [7]. In the latter the contribution due to inhomogeneity to the 
DPCF is obtained from the higher-order direct correlation functions of the isotropic phase. 
The leading term of this expansion gives 



>i,r2;[p]) = j 4°^(ri2,ri3,r23;po)(p(r3) - Po)dv3 (3) 



where p{vs) — Po-\- Pbi^s) with ^^(ra) = po Sg^o A*g 6xp(^G • ra) and 03*^^ is the three-body 
DPCF of the isotropic fluid of density po; Pg are order parameters and the sum is over the 
complete set of R.L.V of the crystal. Eq.(3) satisfies the condition that c'^^^ is zero in the 
isotropic phase and depends on the amplitude and phase factors of density waves which arise 



in a crystal due to breaking of symmetry of fluid. Though the higher order terms in (3), 
particularly for softer spheres, are not negligible and will affect the location of fluid-solid 
transition [17], the relative stability of the two cubic structures can however be understood 
on the basis of (3) as the error caused due to neglect of higher order terms in both structures 
are of similar magnitude. 

Barrat et.al [9] have shown that cf^ can be factorized as cf\ri2,ri3,r23) — 
^(^i2)^(?^i3)^(?^23) and the function t(r) can be determined from the relation 

= t(r) J t{rn)t{r23)drs (4) 
Using above relations we find the following relation for c^*^^ (r) : 
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Ji,{\GT)t{T)Bi,{T,G)Yi^%G)]YUv) (5) 



where Cg is the Clebsch-Gordan coefficient,j/(a;) the spherical Bessel function and 
Bi^{r,G) = {ATifj dkkH{k)ji,{kr) j dr'r'\r')3i,{kr')3i,{Gr') 



(6) 



The crystal symmetry dictates that 1 and li+h are even and for a cubic crystal, m = 0, ±4. 
If we write 6^\y) = Z];m C([^)(r)y/,n,(r), the expression given in the square bracket in (5) is 
the expression for C;^"* (r) . The (r) depends on the order parameter and on the magnitude 
of R.L.V. 

In Figs. 1 and 2 we compare few harmonic coefficients (r) of the DPCF of f.c.c and 
b.c.c structures for n=12 for R.L.V of first and second sets, respectively. Wc note that the 

(C) 

values of Qv„ (r) are far from negligible and differ considerably for the two structures; the 
difference is both in magnitude and in r dependence.lt is this difference that plays crucial 
role in giving relative stability to one crystalline structure over the other. As the magnitude 
of G increases the values of c[^^(r) decreases and after the ninth set of R.L.V values become 
neghgible. 

The functional integration of (1) in density space gives Aea:[p]. In this integration the 
system is taken from some initial density to the final density p(r) along a path in the 
density space; the result is independent of the path of integration [18]. Since the symmetry 



conserving part c*^°^ depends on number density only, the integration is done taking the 
density of the coexisting fluid pi as reference. This leads to 

= AM cir2Ap(ri)Ap(r2)cW(ri, r^) (7) 

where cW(ri, ra) = 2 !^ d\ dX'c^^\r; pi + XX'(po - pi)), Ap(r) = p(r) - pi, AM is the 
excess reduced free energy of the isotropic fluid of density pi and po = pi{l + Ap*) is the 
average density of the ordered phase. Since the functional integrations of c*^^-' have to be 
done in the density space specified by the order parameter pc and the number density po 
we define the path of integration by two parameters A and ^ which vary from to 1. The 
parameter A raises the density from to po as it varies from to 1 whereas parameter ^ 
raises the order parameter from to po as it varies from to 1. This integration gives 

= / dvr J dr2P6(ri)pb(r2)c(^nri, r^) (8) 



where 



(ri,r2) =4 /'da f'de f' dXX /' dA'c^^^ri, ra; AA'po; ^Vg) (9) 
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The free energy functional for a crystal is the sum of Ad, ■ 

The grand thermodynamic potential defined as —W = A — Pp J drp{r), where p is the 
chemical potential, is used to locate transition as it ensures that the pressure and chemical 
potential of two phases remain equal at the transition. The transition point is determined 
by the condition AW — Wi — W — where Wi is the grand thermodynamic potential of the 
fiuid. We calculate the ideal gas part of Al^ using the Gaussian ansatz [19] for the solid 
density p(r) = (a/Tr)^'^^ X^j^ exp(— Q;(r — Ri)^), a being the localization parameter and for 
the excess part the Fourier form with pc = exp{—G'^/Aa). 



AW 



1 - (1 + Ap*)[^ +lnp, - ^ln(a/7r)] - ^Ap*¥°)(0) -W Ipaf^^'KC) 



N ^ -^2 2 - - 2 ^ ^ 2^^, 

-^EE/^Gi/^-g-g.^^''^(Gi + Jg) (10) 

^ G Gi ^ 



where c^^\G) ^ pi J c("\r)e'^-^dr and 5^''^(Gi + iG) = pi Eim IciJ''Kr)e'^'''+'^''>''YUr)dr 
We used the above expression to examine the relative stability of f.c.c and b.c.c structures 
under conditions of fluid-solid coexistence, as determined by MC simulations for soft spheres 
with n=12, 6 and 4. In table 1 we give the values of ideal, symmetry conserving and 



symmetry broken contributions to AW/N. We note that the contribution arising due to 
symmetry broken part of DPCF is far from neghgible and its importance increases with 
the softness of potentiaL While it is about one-fourth of the symmetry conserving part for 
n—12, for n—A it increases to nearly one-half. As the contribution is negative it stabilizes 
the sohd phase. Without it the theory strongly overestimates the stability of fluid phase 
specially for the softer repulsions (n=6 and 4) [11,20]. In agreement with simulation results 
we find that for n=12 the f.c.c structure on freezing is favored while for the softer repulsions 
(n=6 and 4), b.c.c structure is favored. From the values of AW/N given in the table we 
also conclude that while (3) is a good approximation for c^'^^ for n > 12 but it overestimates 
the value of c^^^ as the softness of the repulsion increases. Because of this the crystal phase 
becomes stable at lower values of 7 than found by simulations for n = 6 and 4. 

In conclusion, we developed a free-energy functional for a crystal that contains both 
symmetry conserved and symmetry broken parts of the direct pair correlation function. 
We calculated the symmetry conserving part of the pair correlation functions using the 
Ornestein-Zernike equation and the Roger- Young closure relation and used a perturbation 
expansion for c*^^). The c('')(ri,r2) has been expressed in the Fourier series in the center 
of mass variable with coefficient c^*^^ (r) which is function of difference r(= ri — and is 
found to differ considerably both in magnitude and in dependence on r for the b.c.c and 
f.c.c structures. In agreement with simulation results we find that for n=12 the freezing is 
into closed packed f.c.c structure while for soft repulsions (n = 6 and 4) b.c.c phase becomes 
more stable. The predictive power of the free-energy functional developed here can further 
be improved by improving the accuracy of c*^^^(ri,r2) which can perhaps be achieved by 
including one more term in (3) [17]. 
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TABLE I: Ideal {^Wid/N), symmetry conserving (AWo/iV) and symmetry broken (AWfe/iV) 
contributions to /SW/N for f.c.c and b.c.c phases of three inverse power potentials (n = 12, 6 and 
4) at fluid-sohd coexistance. The coexistance parameters 7s, 7;, and L (the Lindemann parameter) 
are also given 





structure 




AWo/A^ 


AW^b/AT 


AW/N 


n = 12,L = 0.15 


f.c.c 


2.80 


-2.29 


-0.50 


0.01 


7, = 1.19, 7z = 1.15 


b.c.c 


2.88 


-2.24 


-0.55 


0.09 


n = 6,L = 0.17 


f.c.c 


2.38 


-1.80 


-0.62 


-0.04 


7s = 2.33, 7, = 2.30 


b.c.c 


2.46 


-1.87 


-0.69 


-0.10 


ri = 4,L = 0.18 


f.c.c 


2.19 


-1.60 


-0.72 


-0.13 


7, = 5.75, 7, = 5.72 


b.c.c 


2.28 


-1.75 


-0.82 


-0.29 
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